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Abstract. Motivated by the paraxial narrow— angle approximation of the Helmholtz equation in do- 
mains of variable topography that appears as an important application in Underwater Acoustics, we 
analyze a general Schrodinger-type equation posed on two-dimensional variable domains with mixed 
boundary conditions. The resulting initial- and boundary-value problem is transformed into an equiva- 
lent one posed on a rectangular domain and is approximated by fully discrete, L 2 -stable, finite element, 
Crank— Nicolson type schemes. We prove a global elliptic regularity theorem for complex elliptic bound- 
ary value problems with mixed conditions and derive L 2 -error estimates of optimal order. Numerical 
experiments are presented which verify the optimal rate of convergence. 



1.1. The physical problem. The standard narrow-angle Parabolic Equation (PE) in three space di- 
mensions is the following Schr6dinger-type equation 



that models the long-range sound propagation in the sea, and is used in the context of underwater acous- 
tics as the paraxial and far-field approximation of the Helmholtz equation in the presence of cylindrical 
symmetry, cf. [25J [10] . Here, r max > r > r m [ n > is the horizontal distance from a harmonic point 
source placed on the z axis and emitting at a frequency /o- The function -0 = ip{r,z,9) depending on 
range, depth and azimuth measures the acoustic pressure in inhomogeneous, weakly range-dependent 
marine environments. The depth variable z > is increasing downwards while the azimuth varies in the 
interval [0 m i n) 0max]; fco = 2 ^ is a reference wave number, the constant cq is a reference sound speed, 
n{r, z, 6) = C{)/c{r, z, 9) is the refraction index and c(r, z, 0) is the sound speed in the water. The bottom 
topography, being variable, is identified in cylindrical coordinates by a positive surface z = s(r, 9). 
For a fixed range r £ [r mul , r max ], we define the r-dependent space domain: 



where obviously, dQ(r) = U^w^r) for ur(r) := {(0,0) el 2 : <= [0 min , max ]}, w 2 (r) := {(z,0 min ) e 
R 2 : z€ [0, s(r, min )]}, w 3 (r) := {(s(r, 0), 0) G K 2 : e [0 mi „, max ]}, and w 4 (r) := {(«, max ) el 2 : ze 
[O,s(r,0 max )]> (cf. Figure Q}. 



The horizontal sea surface of the naval environment is assumed to be perfectly absorbing, so a free- 
release condition i\) = is imposed on u>i(r). We also set tp — on the minimum and maximum 
azimuthal values i.e. at u>2(r) U u!i(r). We denote by wn(r) := wi(r) U W2(r) U u)i(r) the piecewise 
linear boundary segment where these homogeneous Dirichlet conditions are imposed. The acoustically 
rigid bottom is mathematically modeled by the Neumann boundary condition J^- = along the bottom 
surface z = s(r, •), i.e., the variable boundary segment w 3 (r) of Q(r). Even in one space dimension, 
the well-posedness of the standard narrow-angle Parabolic Equation (jl.ip with Neumann condition was 
proved under the assumption that the bottom topography is strictly monotone, cf. pQ. Considering 
the same problem, in [2l [8] the authors verified numerically that significant instabilities develop even in 
strictly monotone downsloping bottom profiles. 
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1. Introduction 



(1.1) 




n(r) :={(z,0)eR 2 : e [0 min , max ], z € [0, s(r, 0)]}, 
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Figure 1. The range dependent domain f2(r). 



Abrahamsson and Kreiss in [TJ [5] proposed alternatively the use of a Robin-type condition as an 
approximation of the Neumann one that yields a well-posed initial and boundary value problem when the 
domain topography is variable. This approximate condition in two dimensions has the form (cf. [24] ): 



(1.2) 



4>z 



so 



ipe = ikos r ip at z — s(r,9). 



We impose (|1.2I) at z = s, we set in a :— f3^(r, z, 9) := ^(n 2 — 1) and arrive at the following 
initial and boundary value problem (ibvp) of Schrodinger type: 



Vv = i div(D a V?/>) + iA/,-0 
V; = 



in S, 

on u!d(t) Vr € [r- 



min : ' maxj i 



(1.3) 



20: 



mini ' max J 5 



posed on the non-cylindrical domain S '■= U re r, rn 



on Cd3(r) Vr 6 [r n 
on Q(r min ), 

. max ]51(r). Here, the gradient is with respect to the 
:. v variables, D a :— ( a ) , r; s = — is the vector normal to the surface z = s and the initial 

\0 a/r 2 J V 1+S e 

condition ipo models the acoustic source. 

Remark 1.1. In view of the ibvp (|1.3I) . we observe that the same term -D a VV> appears at the equation as 
well as at the left-hand side of the Abrahamsson-Kreiss Robin condition. 

1.2. Change of variables. The focus of our interest herein is to write the problem into an equivalent 
form posed on a cylindrical domain where simpler stable numerical schemes can be applied. This is 
achieved by a horizontal change of variables combined with an exponential transformation. Specifically, 
we let 



(1.4) 



y = z/s(r, 9), v(r, y, 9) = e -^ e ^(r, z, 9), 
n(r) ^ D := (0,1) x (6 

mill) ^max) > 

S c— ^ [r mm ,r max ] X D, 



where q(r, 9) = — | In s(r, 9) [531 [BJ. With this choice of q, the initial- and boundary- value problem (|1.3[) 
takes the following form (see [BJ for the details): 



v r = i div(Z? a Vv) + y— v y + i/3 v v in [r min , r max ] x £>, 



(1.5) 



v = 
t 



mm i ' max] 



77 (D a Vu) = i7bc« 
v(r m in,y,0) = v (y,0) 



at y = V(r, 6>) e [r, 
at ?/ = 1 V(r,0) € [r m i,,. r 

VM)e®, 



'mm; ^maxj ? 



X [0 n 

: ■ y 1 1 1 kJ X [^minj ^max]i 
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where7 ? :=(l,0)*, 7bc (r,0):=i[f +i^(f) 2 ],and5 a := ^ JJ with 

A, :^/3^ + ^ 3sg g Sw — ifg, and Uofj/,0) = V s ( r ™n,0)u (ys(r min ,6),9). 

We note that D a is a real, symmetric and positive definite matrix and therefore det(D a ) > 0, [6]. 
Furthermore, due to the definition of q the coefficient of v y in the first equation is a real function, which 
at y = 1 equals to 2Re7b c - 

Certain three-dimensional effects have been observed to influence the acoustic transmission in variable 
domains mainly because the refraction index depends on r, z, and since significant reflections may 
occur between the bottom and the see surface (cf. [19l [14j [27l [T2l [13] ) . In [24], F. Sturm considered 
the Narrow-angle parabolic equation with the Abrahamsson-Kreiss condition in three dimensions over a 
variable bottom in the case of a multilayered fluid medium. 

The single layer case in the presence of azimuthal symmetry where the physical problem is posed 
on one-dimensional variable domains has been analyzed rigorously in J5[ [5]. More specifically, in [5] 
the authors constructed finite difference schemes and proved optimal rate of convergence. In [5], error 
estimates of optimal order in the L 2 - and -norms have been proved for semidiscrete and fully discrete 
Crank-Nicolson-Galerkin finite element approximations. Discontinuous Galerkin methods for the linear 
Schrodinger equation Dirichlet problem in non-cylindrical domains of R m , m > 1 , were analyzed in [9] . 
When m = 1 the resulting problem is the standard Narrow-angle parabolic approximation modeling an 
acoustically soft bottom; for this case the authors investigated theoretically and numerically the order 
of convergence using finite element spaces of piecewise polynomial functions. The Wide-angle parabolic 
equation consists an alternative approximation model of Helmholtz equation in underwater acoustics; for 
a rigorous numerical analysis and numerical experiments on this model cf. [3j [4] [16] . 

1.3. Generalization: The mathematical problem. Motivated by the properties of the physical prob- 
lem, for the sake of a more general mathematical setting, in our analysis we consider the following initial- 
and boundary-value problem of Schrodinger type with variable coefficients and mixed boundary condi- 
tions (Dirichlet-Robin) 

u r = i div(_DVu) + bWu + i/3u + F in [r m j n , r max ] x 



u = at y = V(r, 0) e [r min , r max ] x [0. 



^' 6 ^ rf(DVu) =i\u at y = 1 V(r, 0) G [r min ,r max ] x [0 min ,0 ma:< 



mm? ^maxj ; 



u(r miD ,y,6)=u (y,6) V(y,0)e3). 

Here, $ = (0,1) x (0 min ,0 max ), r, := (1,0)*, while = p(r t y,9), F = F{r 7 y,6) and A = A(r,0) arc 
complex- valued functions. 

For the rest of this paper, we shall assume that the following conditions are satisfied: 

(1.7) D = D(r, y, 0) is a 2 x 2 real, symmetric matrix with det(D) > Vr, y, 0, 



is real, 



(1.8) b=(b 1 {r,y,6), & 2 (r,y,0)) is 
and 

(1.9) &i(r, 1,0) - 2ReA(r,0) <0 Vr,0. 

Remark 1.2. Since £> e R 2x2 , the condition {12]) gives equivalently that D is either positive or negative 



definite for any r, y, 0, which in turn relates to the ellipticity of the operator div(DV-). 



Remark 1.3. As we shall prove later, the conditions (jl.8j) and (|1.9j) are sufficient for L 2 -stability, while 
when (|1.9[) holds as equality the problem is if ^stable also, cf. Theorem 13.11 and Remark 13.21 



Remark 1.4. The form of the Robin boundary condition, considering only the first order terms, is related 
to the elliptic regularity of elliptic problems with mixed Dirichlet-Robin conditions in two dimensions 
proved in Theorem l4.3l The autonomous Section 4 of this paper presents a detailed proof of this argument. 
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Remark 1.5. The acoustic problem (|1.5[) is a specific case of the problem (|1.6j) for u := v, uq := Do, 
D := £> a , b := (y^, 0), (3 := p v , F := and A := 7 bc , satisfying (|L~7|) . ([TS]) . and JOJ) as equality, [6]. 

1.4. Main results. The problem analyzed here is motivated by an important physical application. Nev- 
ertheless, the general mathematical setting encompasses the very interesting aspect of approximating 
numerically a multi-dimensional ibvp of Schrodinger-type with mixed conditions and coefficients depend- 
ing on the evolutionary variable. 

In this paper, we apply the Galerkin method on the general problem (|1.6p using piecewise polynomial 
finite element spaces. We construct fully discrete Crank-Nicolson-type schemes in r for which we prove 
stability and optimal rate of accuracy in the L 2 -norm. Numerical verification of the optimal rate of 
convergence is also presented. 

The weak formulation of the problem is presented in Section 2. We define an appropriate r-dependent 
sesquilinear form which is, in general, non-Hermitian. As it is common, the rate of accuracy is investigated 
by using certain properties of the projection induced by this form. The projection being r-dependent and 
the fact that a two-dimensional r-dependent Robin boundary condition appears in (jl.6[) make the analysis 
difficult. We estimate the projection error and its r-derivative in the H 1 - and L 2 -norms (cf. paragraph 
2.3, Propositions l2.3H2.5p . The later is accomplished by applying an Elliptic Regularity Theorem for 
two-dimensional complex boundary value problems with mixed Dirichlct and Robin conditions, proved 
in Section 5. In the proof of Proposition 12.51 where the r-derivative of the projection error is estimated 
in the L -norm, we present a very refined argument when treating the boundary terms. 

In Section 3, we write (|1.6p in a weak form and prove L 2 -stability, and i/ 1 -stability in the case where 
(|1.9p holds as equality, so that the sesquilinear form is Hcrmitian. We then construct a fully discrete 
Crank-Nicolson scheme in range r that is shown to be L 2 -stable. Even though the evolutionary variable 
is discretized by a standard Crank-Nicolson method, the error analysis presented in this section is non- 
standard. This is due mainly to the fact that the form and the projection used are r-dependent and 
calculated at the mid-points of a uniform range partition. We define properly a test function split in 
two terms involving projections applied on second order derivatives (cf. Remark 13. 6p . use the projection 
estimates of Section 2, and derive an optimal error estimate in the L 2 -norm. 

A general complex elliptic boundary value problem posed on a two-dimensional rectangular domain 
with mixed boundary conditions is analyzed in Section 4. If Dirichlet or Neumann conditions hold along 
the boundary, then in the weak formulation of the boundary value problem the trace integral terms 
vanish. A general approach of proving global regularity, [18], is to prove this estimate for half-balls, 
and then by change of variables, stretch the compact boundary locally and cover it by a finite union of 
half-balls. In our case, we analyze a complex elliptic problem posed on a rectangular domain of M 2 . The 
boundary is compact and consists of four linear segments along which Dirichlet and Robin conditions are 
imposed. We apply directly on this domain the half-balls technique without change of variables as the 
boundary is already stretched locally. Further, we define appropriate test functions, in order to eliminate 
the trace terms from the weak formulation of the problem and prove the regularity estimate in Theorem 
14.11 The result is extended in Theorem 14.31 Our proof covers a class of Robin conditions related to the 
coefficients of the pdc of the boundary value problem, a special case of which is the Abrahamsson-Kreiss 
condition of underwater acoustics. 

Finally, in Section 5 we report on the results of some numerical experiments performed with our 
method, verifying experimentally the optimal order of convergence. 

2. An elliptic projection 
2.1. Preliminaries. Let 2) = (0, 1) x (#i,#2)- For r in [r m i n ,r max ] fixed, we define 

dy=vo ■= {(r,Vo,0) e K 3 : 9 € [6i,6 2 ]}, d e =e := {(r,y,0 o ) G K 3 : y € [0,1]}, 

and denote by 7? 1 (£>) the associated usual (complex) Sobolev space. In order to deal with the Dirichlet 
boundary condition we shall make use of the space 

H%(P) = {u e H 1 ^) : u\ 9Sd = 0}, 
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where 8Dd ■= (d y =o) U (80=0^ U (de=g 2 ). Hq(D) is the space of functions in if 1 (33) which vanish on y = 
0,6 = 6i,6 = 9 2 . We denote the L 2 (33) inner product by (•, •) : L 2 (33) x L 2 (33) -> C. ||tt|| := ( J s |tt| : 



denotes the induced L -norm, while Hulli := [\\u\\ 2 + \\u y \\ 2 + \\ug\\ 2 j is the usual H 1 ^) norm. 

Let d®n :— d y —i be the part of <933 where the Robin boundary condition of problem (jl.6l) is posed, 
and let 

< u,v >:= I u(l,6)v(l,6)d6, 



denote the inner product on H 2 (<933fl). In addition, we shall make use of the norms: 

2 v£Hl(®):v\ov R =g 
II \<U,V>\ 

\ u \-h,dT> R ■= SU P ii^Ti • 

veH 1 ^),^ \\ v \\l 

Let t <E N and 5/, be a finite dimensional subspace of Hq(D) consisting of complex- valued functions 
that are polynomials of degree less than or equal to r in each interval of a non-uniform partition of 33 
with maximum length h £ (0, h*]. It is well-known, that the following approximation property holds: 

inf {\\v-x\\+h\\v- X \\i} <Ch s+1 \\v\\ s+1 , Vv£H s+1 (T>), 

(2.1) * eSh 

Vft,e(0,/i*], s = o,...,r. 

Also, we assume that the following inverse inequality holds: 

(2.2) < C h' 1 \\(f,\\ V<P£S h , Vh£(0,K], 
which is true when, for example, the partition of 33 is quasi- uniform, 

2.2. Definition of a sesquilinear form. Without loss of generality we assume that D is positive 
definite. For any r in [r min , r max ] we define the sesquilinear form B(r; •, •) : Hq (33) x (33) — > C 

B{r;v,w) :=(D(r)V«, Vto) -if / A(r, 0)u(l, 0)to(l, 0)d0 - (6(r)Vu, to)) 
+ S(v,w), 

for 5 a sufficiently large positive constant. Obviously, it holds that 

(2.4) |B(r;«,ic)|<c|M|iHli, 

for any v, to £ Hq(D), uniformly in r. 

We observe that (DVv, Vv) > co||Vu|| 2 for a constant Co > uniformly in r and v, since D is real 
symmetric and positive definite. Therefore, by the trace inequality we obtain for v £ Hq(D) 

ReB(r; v, v) ={D{r)Vv, Vv) + S\\v\\ 2 

+ Im( J \(r, 6)v(l, 6)v(l,6)d6 - (b(r)Vv, w)) 

>c ||Vt;|| 2 + 6\\v\\ 2 -c?|H|||«||i. 
Thus, by choosing 5 sufficiently large, it follows that there exists positive constant C such that 

(2.5) ReB{r;v,v)>C\\v\\l 
uniformly, for any r and any v £ Hq(D). 

Remark 2.1. If D is negative definite we may use 



B n (r; v, w) := - (D(r)Vv, Vio) X(r, 6)v(l,6)w(l, 9)d9 - (b(r)Vv, to)) 

+ S(v,w), 
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with S a sufficiently large positive constant. In this case (|2.4p and (|2.5[) also hold since now —D is positive 
definite. 

2.3. Projection estimates. Let Rh{r) : Hq(1)) — > be a projection operator defined by 

(2.6) B(r; R h (r)v, <j>) = B(r; v, cf>) e S h . 

Obviously, since (I2.4[) and (I2.5[) hold true, then by Lax-Milgram Theorem the projection is well defined. 
Let us now define the operator 

C*(r)w := -div(DVw) + ibS7w + i[b ly + b 2 e]w + Sw in D, 

(2.7) w — on 9Dd, 
if{DWw) = i\*w on &S R , 

with A* a complex- valued function to be chosen appropriately in the sequel. For (f> € H^T)) we get 



(£*(r)w,<j>) = (L>Vw,V0) - / i\*w(l)<j>(l)de 

+ i(bVw, 4>) + \{[b ly + b 28 ]w, <p) + S(w, 4>) 

= (DVw,V4>)- f i\*w(l)0(l)d8 



(2.8) " m,n 

- i([h v + he]w, 0) - i(bw, V0) + i / 6i(l)u;(l)0(l)d0 



+ i([h y + b 2e ]w, 4>) + 5(w, 4>) 
= {D\7 W , V</>) + i / [&i(l) - \*]w(l)4>(l)dd - i(bw, V0) + 5{w, < 
Since I?, b and 5 are real, then for any <f> in Hq(D) it follows that 



(£*(r>, </>) = (DV<P, Vw) - i J [h(l) - \*}<t>(l)w(l)d6 
+ i(&V</>, w) + 8(<j>, w). 

Setting 

(2.10) A*:=6i(l)-A, 
we obtain 



(2.11) {L*{r)w,4>)=B{r;<l>,w), 
and thus 



(2.12) (£*(r)w,®=13(r;ct>,w), 

for any <f> £ Hq(T>). Throughout the rest of this paper, we consider A* given by p.!0|) . 

Remark 2.2. We observe that in the case of the specific problem (|1.5[) . b\ — y^-, b 2 — 0, A = | 

— «- 1 [ s r . a ( se\ 2 ~\ 

A. 



s 2 L s r z \ s / 

Proposition 2.3. There exists a positive constant c such that if v 6 ffg(5D) n iP(!D) £/ien 

(2.13) ||i2fc(r)i;-«||i < c/i T ||u|| T+ i, 
and 

(2.14) ||i2 h (r)« - w|| <c/ i r+1 ||«|| r+1 . 
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Proof. We set e := Rh(r)v - v, use (|23j) . (|2~4| and (f27T]) to obtain for e S h 

c\\e\\l < ReB(r; e, e) = Re£(r; e, fl/,(r)u - v) = ReB{r; e,<p-v) 

< cHelK inf ||^-«|| 1 <c||e|| 1 fe'-||t;|| r+1) 

4>es h 



which establishes ()2.13|) . 

Let now w be the solution of the problem: C*(r)w — e. Then by using (|2.12j) . the approximation 
property and elliptic regularity, proved in Theorem 14.31 we get for <p G Sh'- 



||e|| = (C*(r)w, e) = B(r; e, w) — B(r; e, w - 



< cllelli inf lilt; 



< ch T \\v\\ T+1 h\\w\\ 2 < ch T+i \\v\\ T+1 \\el 



which yields (|2.14p . 



□ 



Proposition 2.4. Let r e C 1 ([r m i n , r max ], H r (®)). Then it holds that 

(2.15) \\d r (R h (r)v(r) - vfr))^ < C h T ( ||«|| T+1 + ||SHIr+i ) . 

Proof. We set e :— Rh{r)v(r) — v(r). Let v : [r m i n ,r max ] — ► H r {T)) and e(r) = Rh(r)v(r) — u(r) for 
r € [r min ,r max ]. Then, we have 

B(r;e(r),0) =0 V^e5,, 
Differentiating the above relation with respect to r we obtain 

B(r; e(r), 0) + £(r; e(r), 0) = e S h . 

Now, for cf> E Sh we have 

c\\e{r)\\\ < B{r; e(r), e(r)) = B{r; e(r), e(r) + </>)- B(r; e(r), 0) 
= B(r; e(r), e(r) + </>) + B(r; e(r), </>) 

<c[||e(r)||i||e(r) + 0||i + ||e(r)||i||0||i 

< c[||e(r)||i||e(r) + 0||i + ||e(r)||i(||e(r) + <f>\\i + \\e(r)\\x) 

"(11^)11! + ||e(r)|| 1 )||e(r-)+0|| 1 + ||e(r)|| 1 ||e(r-)|| 1 



= c 

< c 

< c 



(\\e(r)\h + ||e(r)||i) inf \\d r (R h v)(r) - 8^ + ^ + Mr^Wetfh 

4>es h 

(l|e(r)|| a + \\e(r)\U) inf - + lleWIKHeWHx 

0eSh 



c||e(r)||i ||e(r)||i + inf ||<9 r u- 

L <pGS h 



+ c||e(r)||i inf \\d r v ■ 
<pes h 



The claim of the proposition follows by using the approximation property 12.11 with s = r + 1. 



□ 



Using a technique introduced in |17j . we are able to show the following optimal order approximation 
result for the time-derivative of the elliptic projection. 

Proposition 2.5. There exists a positive constant c such that 

(2.16) ||9 r (R h (r)v(r) - v(r))\\ < C h T+1 ( ||«|| T+1 + \\d r v\\ T+1 ) . 
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Proof. We set e := Rh(r)v(r) — v(r). Let w be the solution of the problem: C*w = e. For \ <= &h we 
have 



\e(r)\\ 2 = (C*w,e(r)) = B(r; e(r),w) 



= Re 
< c 
< 



B(r;e(r),w- \) - B{r;e,x) 

- xlli + l|e(r)||i||u> - xlli 



- Re 



B(r; e, w) 



c(p(r)|| 1 + || e (r)||) inf \\w - xlli - Re\B(r;e,w 
< c/i r ^||w|| r+ i + ||9 r w|| r+ i^/i||w||2 — Re B(r; e, w) 



For convenience we set I := Re 



B(r: e, w) 



. First, observe that 



a $m ax v 
d r X(r, 8)e(r, 1, (9)uJ(r, 1, 6>)d0 - (<9 r 6Ve, w)J , 



so that 



7 = Re 



(dr-DVe, Vra) - i 



<9 r A(r, 0)e(r, 1, 6>)uJ(r, 1, 0)d0 - (d r bVe, w) 



By the definition of the inner product < u, v > we have 



h := Re 



i(<9 r 6Ve, w) 



= Re 



- i(9 r [6i s + b 2 e]e, w) — i(d r be, Vro) 



+ i 



(fl^.6i)(l)e(l)«J(l)dfl 



< c||e||||u;||i + Re 



i < d r bie, w > 



We set 7 2 := Rc 



i / e max d r \ewd9 . Using the estimates above we obtain 



h+h< Re 



i < [d r \ — d r bi]e, w > 



+ c \\ e \\\\ w i 



< c|e|_i OT JM|i + c||e||H|i < c 



to i- 



In addition, 



J 3 := Re 



= Re 



j]{d r DVwt)de - (div(d r DVw), e) 



< c|H| 2 ||e|| + Re 



< rj(d r DVw), e > 



<c|Ml2l|e|| + |e| i 



\w\\ 2 < c 



e -f e-i 



W 2, 



so that 
(2.17) 



Kc 



e + el i 



,8X>b 



to 



Now, for g e iJ 2 (dDii) we consider the elliptic problem 

- div(UVz) + ibVz + i[b ly + b 2 e]z + Sz = in D, 
z = on 

if(DVz) = \\*z + g on dD R . 



Then we have = B(r; e(r),z)— < g,e > and thus, 



< g, e >= B(r; e(r), 2) = B(r; e(r), 2 - 0) V0 e S h . 
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It follows then that 

I < e,g > | < c||e||i inf ||z-<£||i, 
4>£S h 

and therefore, 

| < e,g > | < c/i T ||w|| T+ i/i||z|| 2 . 
The elliptic regularity result (cf. Theorem 14.31 and Remark 14.4)) for the solution z of the elliptic problem 
above, reads 

NI2 < cb|i,as R - 

Thus we have 



el 1 



sup 



< e, v > 



< ch T+i \\v\\ T 



+17 



+1 > 



ueH 1 ^)^ \\ v \\i 

and subsequently, using the elliptic regularity of w, cf. again Theorem 14.31 we arrive at 
||e(r)|| 2 < ch T+1 (\\v\\ T+1 + ||a r «|| T+1 )||w|| 2 + c[||e|| + \e\_ 1/2tdVR \\\w\\ 2 

< c|| W || 2 /i T+1 (||i;|| T+1 + ||flHlr+i) < c||e(r)||/i T+1 (|M| T+ i + \\d r v\\ T 
which completes the proof of the proposition. 

3. A Crank-Nikolson-type FULLY DISCRETE scheme 

3.1. Weak Formulation. Let £ Hq(D). Multiplying the partial differential equation of (11.61) by 

and integrating by parts we have 

(u r (r),<f)j =i - (.D(r)Vu(r), V<pj + j rfD{r)Vu{r)4>ds + (o(r)Vu(r)> 



□ 



+ i(/3(r)«(r),0) +(F(r),0) 



(3.1) 



D(r)Vu{r),V(t)j -i 
+ i[/3(r)u(r),A +(F(r),4»: 



A(r, 9)u(r, 1, 9)<p(l, 9)d9 - (b(r)Vu(r), } 



= -is(r; u(r), </>) + i((/3(r) + S)u(r), <t>) + (F(r), 0), 
for any </> G ZIq(ID). In the following theorem we prove that (|3.1|) defines u in iJ^S)) uniquely. 
Theorem 3.1. The weak problem (|3.1|l /tas ai most one solution in Hq(D). 

Proof. Let m e -ffo(®) be a solution of (|3.1|) . We set </> = u in (|3.1|> . integrate by parts, use the facts that 
D is a real, symmetric matrix, that b is real, and take real parts. More specifically, we obtain first 



(3.2) 



Observe that 



(u r , u)=-i (DVu, Vu) - i( / A(r, 9)u(r, 1, 9)u(l, 9)d9 - (bVu, u)) 
+ i(0U,u) + (F,u). 



Re(b\7u,u) = -^{b ly u,u)~^{b 2e u 7 u) + ^ / h(r, 1, 9)\u(r, 1, 9)\ 2 d9, 



since is real and u = at y = 
arrive at 

1 d 



2dr 



u\\ 2 = 



ii^max- Since I? is real, then using this observation in (|3.2[) we 



[-ReA(r, 0) + -b x (r, 1, 0)]|u(r, 1, 9)\ 2 d9 



-([hy + b 2 e]u, u) ~ (Im(/3)u, u) + Re(F, u). 
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Using the condition (|1.9|) and Gronwall's inequality we obtain the stability estimate 
(3.3) \\u\\<c\\u \\+cJ^\\F\\dr. 

Uniqueness of the solution u follows readily from the estimate above. □ 



Remark 3.2. If f| 1 . 9[) holds as equality then the sesquilinear form is Hermitian. Therefore, if F = 0, using 
()3.10 . setting <f> = Ur — + S)u and taking imaginary parts we obtain 

-— ReB(r; u(r), u(r)) < ci||u||i < cReB(r;u(r),u(r)), 
so that c||u||f < KeB(r;u(r),u(r)) < c||uo||f, i.e. we also obtain an H 1 stability estimate. 

Remark 3.3. Note that for the specific case of problem 1)1.50 . if A/> is real, we have F = 0, Im(/3) := 
Im(/3„) = — |j, b\ — y^-, b 2 — and 1)1.90 holds as equality, therefore (cf. the proof of the previous 
theorem) we obtain the conservation property 

\M = \\vo\\ 

for any r, while the problem is also if -stable. 

3.2. The numerical scheme. For N > 1 integer, we consider a uniform partition in range r m i n = r° < 
r 1 < ■ ■ ■ < r N = r max , < n < N, k := r n+1 ~ r n = ± for any n < N - 1, and set r n+ i := r "+f +1 . If 
u is the solution of the continuous problem (11.60 . we approximate u(r n+1 ) by U n+ G as follows: for 
J7" known we seek U n+1 G Sh such that 

-U n+1 -U n ,\ U n+1 + U n ,\ . „[/" +1 + C/ r ' 

(3.4) 

+ (F(r"+3),0), 



for any ^ G Sh, and any < n < N — 1. In order to obtain an optimal order approximation we shall take 
U° := R h (r°)u eS h . 

Remark 3.4. Let w G Hq(D). The need of the condition 1)1.90 appears once again (recall that (11.90 was 
used for the L? stability of the continuous problem). More specifically, since w G Hq(D) then by 1)1.90 
the following inequality holds 

(3.5) Re{ -m(r;w,w)\ < c\\w\\ 2 . 

Theorem 3.5. The fully discrete scheme {3.$ is L 2 -stable. 

Proof. In p. 40 we set 4> — U n+1 + U n G Sh C Hq(T>), take real parts and use the estimate of Remark 
13.41 to arrive at 

(1 - ck)\\U n+1 \\ < {1 + ck)\\U n \\ + ck\\F(r n +i)\\. 
Choosing k sufficiently small we get a stability result for the scheme p. 40 

\\U n \\ < c\\U°\\ + cmax||.F||, 

n< N 

for any 1 < n < N. Consequently, uniqueness of solution in Sh is also established. □ 
3.3. Error estimates for the fully discrete scheme. 

3.3.1. Preliminaries. We define in Sh C i?o(2)) t ne quantities 

i k 2 i i 

ffl := U n - R h (r n+ i)u(r n ) + —R h ( r n+ 2)u rr (r n+ ^), < n < N - 1, 

8 

gg+l un+ i _ Rh ( r n+^ u ( r n+l^ + & ^(r^+S ) Urr ( r »+f ), < Tl < JV - 1, 

8 

where J7™ is the solution of the fully discrete scheme 1)3.41) . 
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Remark 3.6. The main idea is to mimic the continuous problem. In the fully discrete scheme, we set 
4> := #2 +1 + ®\ as test function. The choice of $2 +1 1S n °t standard and is made in order to treat 
efficiently the r-dependent sesquilinear form at the midpoints of the partition and since the projection is 
range-dependent. Therefore, in </>, the projections are computed in r" + 2 . The introduction of the specific 
additive term 



(3-6) 

is motivated by the approximation 



k' 2 



R h {r n+ 2)u rr (r n +-2) 



u(r n ) + u(r n+1 ) ,i fc 2 ,i 4 

2 --u(r 2 ) = —u rr {r n+ 2) + 0{k i ), 

used in (|3.10[) . The residual being of order 0(fc 4 ) permits us to apply then in p. lip the inverse inequality 
without loss of optimality in space, and avoid thus any integration by parts (denote that in this case 
suboptimal trace integral terms would appear, as the problem is posed in Sh C Hq(D) ^ Hq(D)). 
Furthermore, the term (13.61) is related to the approximations 

u(r n+1 ) - ^-u rr (r n+ i) = u(r n+ ^) + ^U r (r n+ ^) + 0(k 3 ), 
8 2 

u(r n ) - —u rr {r n -^) = u{r n -^) + -u^r* - *) + 0(k 3 ) 
8 2 



used in the proof of Lemma 13.81 when treating the r-derivative of the projection error. 
We notice that 

jjn+X _ jjn ^gn+i _ q„ R h (r n +h) u (r n+1 ) - R h (r n+ ^)u{r n ) 



k k 
jjn+i + jjn gn+i + Q n R^H ) u (r n+1 ) + R h (r n +^)u(r n ) 



Replacing these identities in the fully discrete scheme we obtain 

' R h (r n+ i)u(r n+1 ) - R h {r n+ ^)u{r n ) 



(3.7) 



9^ +1 - eg 
1 

iB[r 



9£ +1 + 0? 



+ i((/3(r n +i) + 6) 
i((/3(r 



R h (r n+ i)u(r n+1 ) + R h (r n+ i)u(r n ) 
2 ' 

^i?,(r"+^) Urr .(r"+^),0) 



9™ +1 + e r i 



.fc 2 
l- 



i, n R h (r n+ hu(r n+1 ) + R h (r n +^)u(r n ) 

•a) +5) '- i L ^— L : 



— ( (p (r n + 3 ) + 6) R h (r n + * )u rr (r n + § ) , <fy + (F(r n + 1), 4 
From the continuous problem we have that 

+ i((/3(r"+5) + 5)u(r n+ ?), + (F(r n+ ?), 



(3.8) 
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We now solve Q3.8P for (F(r n+ a ), 0), replace in (|3.7I) . and use the definition of the elliptic projection Rh 
to arrive at 

'R h (r n+ i)u(r n+1 ) - R h (r n+ i)u(r n ) 



Wr(r™ + 5), 



i ^ W n+i w(r n+1 )+w(r") , n+i , a 



(3.9) 



' l+ *;^i? A (r"+'K r (r™+i)^ 



9£ +1 + ^ 



-i6(r"+2 
+ i£(V 

+ i((/3(r n+ 5)+5) 

- ( (/3 (r n + * ) + J) i? h * ) Urr (r n+ * ),</>) . 

By Taylor's formula the following identity holds for ri,f2 € (r™,r™ +1 ) 
M (r") + M (r™ +1 ) , „, i. k 2 . n+ i. A: 4 

- U(r 2 ) = — Wrr(r )+ 2 .16 4j Krrr(n) + ^rrrr(^2j 



fl/ l (r"+5) u (r"+ 1 ) + i2 h (r"+i)n(r") , 
u{r 2 J 



Using the above in (|3.9|) we obtain 

r^ 1 

iB 



(3.10) 



iB(r n +i;~ [R h (r n +i)u rr (r n +i) - u rr (r n+ ?)} , 



iB 







(r" + 2 ; - — — [u rrrr {ri) + u rrrr (r 2 )},(t^ 
-iB(r n+ i: 



2-16-4! 

Therefore, applying an inverse inequality we obtain 



Urrrr(ri) + U rrrr (r 2 )] 1 



(3.11) 



Re 



„_i_i M(r" +1 ) + tt(r") , „,i, k 2 i i 

i6(r" + 5;-i ; v ; -u(r"+a) fl/ l (r n+ 2)u rr (r" + = ), 1 



< c/c 4 max||u rrrr ||i||(/<||i < c/c 4 /i 1 max||u rrrr ||i||^)|| 

r r 

In addition, the Taylor formula gives 

■R h (r n+ i)u(r n+1 ) + R h (r n+ i)u{r n ) 



= i((/3(r"+5) + 6) i? /l (r"+5) u (r"+5) - u(r n+ ^) 
+ i((/3(r"+3) + «S)fl h (r"+T) Urr (r" + T) + 



w(r"+2) 



i( (/3(r™+5) + S) R h {r n+ ?)u(r n+ ?) - u{r n+ ^) 



2-16-4! 



[■u rrrT .(ri) + u rrrr (r 2 )] 



Thus, we obtain 



(3.12) 



+ i((/3(r" + ^ + 5)(i? h (r"+^)-7) 

/ 1 [k 2 i A; 4 i \ 

+ i[{/3{r n+ 2 ) + S) [yUrr(r™ +5 ) + — — 4 , Krrr(ri) + U rrrr (r 2 )]J , . 

i? h (r"+3) u (r"+ 1 ) + R h (r n+ i)u(r n ) n+i 
u(r 2 J 



Re 

< 



c{h T+1 + k 2 )\\4>\ 
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Also Taylor gives for r 3l r 4 € (r n ,r n+1 ) 



u(r n+1 ) - u{r n ) 



8-3!' 



"(73) + u rrr (r 4 )] + u r (r n+ ^) 



therefore, 



R h (r n +^)u{r n+1 ) - R h (r n+ ^)u(r n ) 



= -{[R h (r n+ i) - I]u r (r n+1 ),^j - ^([R h (r n+ ^) - I]{u rrr (r 3 ) + u rrr (r 4 )}, cj>) 
k 2 / \ 

- ^— -^[urrrir-i) + U rrr (V4 ) , 4>j ■ 



The above yields 
(3-13) Re 



R h (r n+ i)u{r n+1 ) - R h {r n +?)u(r n ) 



u r (r n+ '),^)] <c{^ r+1 + fc 2 }||0||. 



Let us now assume that k < 1 and k < ch? . In (|3.9[) , we take real parts and use (|3 . 11 1) , (|3.12l) , and 
(|3.13[) to obtain 



Re 



(3.14) 



an \ i r 

( ° 2 k 01 A)] <c{h^ + k 2 }u\\ 



+ Re 



-i8r" + i 



L,0)+i(O9(^+*) + 5)^^1,0) 



In the above estimate, we set 
obtain for A; sufficiently small 

(3.15) 



?2 +1 + 0" G Sh c ffo(S), and use the estimate of Remark 



to 



r 2 +1 \ 



< 



l + ck 
1-ck 



■A, < n < N - 1, 



wnere ^ < cfc ^ +1 + fc2 > . 
Let us now define 



(3.16) 
and 

(3.17) 



T := U n - R h (r n )u(r n ), 0<n<N. 



B ( 2 n+1} := -R h (r n+1 )u(r n+1 ) + R h {r n+ ^)u{r n+1 ) - —R h (r n+ ^)u rr (r n+ ^), < n < N - 1, 

8 



B[ n) := ~R h (r n )u(r n ) + R h {r n+ ^)u{r n ) ~ —R h {r n+ ^)u rr {r n+ ^), < n < N - 1. 

8 

So, we obtain 



(3.18) 



9 n+l = gn+l _ B (n+1) Q < T» < iV - 1, 



9 « = 



B| n) , 0<n<JV-l. 



We replace in (|3.15[) so that for any 1 < n < N — lwe arrive at 



9" +1 — £?2" +1 ^ 



1^)1 

1 + ck 



(3.19) 



- £ 



(")| 



< 



1 - ck/ 
1 + ck^ 



- B. 



(")| 



.4 



4 
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3.3.2. The estimates. We prove first the following lemmas. 
Lemma 3.7. For any < n < N — 1 it holds that 
(3.20) 



<ckh T+ \ 



Proof. By the definition of B 2 n+1 \ B± we obtain 



B (n+i) _ B (n) = [ i?/i(r «+i ) _ T\[u(r n+1 ) - u{r n )\ - / d r [R h {r)u{r) - u{r)]dr. 



Using Taylor's theorem, we obtain for r$,rQ <E (r n ,r 



fc 3 



w(r n+1 ) - u{r n ) = ku r (r n+ 2) + ^-^[u rrr (r 5 ) - u rrr (r 6 )], 
The result now follows from the estimates of Rh(r)v(s) — v(s) and d r (Rh(r)v(r) — v(r)). 
Lemma 3.8. For any 1 < n < N — 1 it holds that 
(3.21) \\B { 2 n+1) - B ( 2 n) \\ < ck[h T+1 + fc 2 }. 

Proof. By the definition of B 2 n+1 \ we have that 



(3.22) 



B (n) _ B (n+1) = Rh{r n + ^\k r ^ nH) _ ^ n+1) 



Rh{r n —i) —u rr {r n ~i)-u{r n ) + R h (r n+i )u(r n+i ) - R h (r n )u(r n ) 
l a 



Using Taylor's theorem we have that for r-j € (r n+ 2 , r n+ ), rg € (r™ 2 ; r ra ) 

u ( r «+i) _ ^ Urr (r"+^) = u(r n +i) + ^u r (r n +i) + J^u rrr (r 7 ), 

fc 2 1 1 k 1 fc 3 

w(r n ) - y« rr (r"^) = u(r"-2) + _ Ur ( r "-2) + __ Urrr ( r8 ). 

Replacing these expansions in (|3.22[) it follows that 



R (n) R (n+1) 



= -R h {r n+1 i)\u{r n+1 i) + - Ur (r n+ ?) 
k 



(3.23) 



+ R h (r n -i) u (r"-t) + ^ Ur (r n -2)j + R h (r n+1 )u(r n+1 ) - R h (r n )u(r n ) + B x 



[d r R h {r)u(r) - u r (r)]dr 
1 (d r R h (r){u(r) + -u r {r)] - [u r (r) + -u rr (r)]jdr 



+ u(r n+L ) - u(r n ) - u{r n+ ^) + u{r n -?) - '^u r {r n+ i) + '^u r {r n -?) + Bi, 
where \Bi\ < ck 3 for h < 1. Expanding in Taylor series around r", r n+1 we finally have that 

\u{r n+1 ) - u{r n ) - u(r n+ ^) + u(r n -^) - ^ti r (r n+ 2) + ^r(r n_ 3)| < ck 3 , 
and the result follows from (|3 . 23[) . 

Remark 3.9. Obviously we assumed n > 1, since we used the nodal point r n_ 2 . 



fc 



fc 
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Lemma 3.10. We have 

(3.24) \\B[ 0) \\ <ckh T+1 + ck 2 

(3.25) 

(3.26) \\B { 2 ' 1+1) \\ <ch T+1 +ck 2 , 0<n<N-l. 

Proof. We use the fact that for rg e (r°, r 1 ) 



IB^H < ckh T+1 +ck 2 , 



k 2 



fc 3 



U (r ) = u(r 3 ) - -Mr^ 2 ) + -5-Wrr(»" 2 ) + 3— ^7 «rrr (j"9 ) , 



S 



8-3! 



and obtain 



B 



(0) 



[9,.i?h(r)u(r) - u,.(r)]rfr - - [R h {r 2 )u r (r 2 ) - u r (p)] 



1 k 1 

+ u(r5) - u(r°) - -u r (r2) + B 2 



for |i5 2 1 < cfc 2 . Finally, using 



|u(r2 ) - u(r°) - -u r (ri )\ < ck 2 , 



we arrive at (|3 . 24[) . 

By Lemma 13.71 applied for n = 0, we get 

\\BP - B[ 0) \\ <ckh T+1 , 

so using (|3.24[) the estimate (|3.25[) follows. 
Using that 



\B 



< \\B. 



(n+i) _ B (n) 



\B^ n) -4 n_1) l 



LB. 



Ml 



and Lemma [3^81 we obtain 



\B 



< IIS 



(!) 



ch T+1 + ck z , 



and by (j3\~2"5]) we arrive at (|3T2^) . □ 

We now estimate 1 . 
Lemma 3.11. If k = 0(h) then 

(3.27) H^H <ch T+1 + ck 2 . 

Proof. We use the continuous problem and the fact that 6° — 0, set cf> — 8 1 in the fully discrete scheme, 
take real parts and use the inverse inequality to obtain 

ln +i_ n n x < R h (r n+1 )u{r n+1 ) - R h (r n )u(r n ) 



LB 



(■ 



-1/2 



6> n+1 



^)+i([/j(r"+ 1 /2 )+(5] ^ 



^(^" +1 )^(r" +1 )+^(r").(r") _ u(r „+l/ 2)j ^ 
+ 1(^1/2) + ^^(^Mr^+^^M^) _ u(r „ +1/2)) ^ 
Obviously, if £> has smooth coefficients and g, a, b are smooth, it follows that 

B(r n+1 / 2 ;a,b)\ < ^|s(r n+1 ; a, b) + B(r n ; a, b)\ + cfe 2 1| « || x || fo)) x , 

S(r" +1 ;a,6)-S(r";a,6) < cfc||a||i||6||i, 
,g(r n+1 ) +g(r n ) 



B(r; 



-,6)-S(r; 5 (r"+ 1 /2) ;&) < ck 2 ]^. 
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So, we get 

B ( r n+i/2. R h {r^)u(r^) + R h (r-)u{^) _ <r „ + i/ 2)) 
g^n+l/2. ^(r" +1 )u(r"+ 1 ) + i? A (r") M (r") _ M (r"+ 1 ) + u(r") ^ ^ + ^ 



_ x flfc(r n+1 )u(r" +1 ) + R h {r n )u{r n ) u{r n+1 ) + u(r n ) 



l n ( n R h (r n )u{r^) + fl„(r"Mr") u^ 1 ) + M (r") \ 
jS^r ; ,(pj +A2+A1 

l B (r n+1 - Rh{rn)<rn) U{rn \&\ I 1 B ( r n. R h{r n+1 )u(r n+1 ) u(r^) 

+ A 2 +Ai = 

2 s r ' — 2 2-^)+^ 



+ o B V ' 2 2 +-^4+^2+ A = 

^3 + ^4 + ^2+^1, 



where 



\A 2 \ < cfc 2 ||0||i (since || R h (r)u(r) ||i < c/i T + c), 
|*4 3 |, \A A \ < ck^UlU (since ||i2fc(r)«(r) - t*(r)|| x < ch T ). 
Therefore, we obtain setting n — 0, (j) = 8 1 and using the inverse inequality 

H6* 1 !! 2 < ck[h T+1 + kh T + /c 2 ] 1 1 6* 1 1 1 i + ckp 1 ^ 

< ckhr l [h T+1 + kh T + fc 2 ]!^ 1 !! + ckp^i 2 . 

So for 0(k) = O(h) the result follows. □ 

Lemma 3.12. If O(k) = <D{h) then for any n > 

(3.28) \\e n+1 \\ <ch T+1 + ck 2 . 

Proof. Since k < cbh then the inequality (|3.19p holds. So, using (|3.19p and Lemmas 13. 7H3. 81 we arrive at 

|| r +l _ B («+D|| < c ||0» _ B (n)|| + cA .^+l + cA .3 ; 

therefore, setting £ n := 9" — for n > 1, we obtain 

||£„+1|| < / l+£*\ ||£„|| + + c/fc 3 

V 1 — ckJ 

< (i±^) 2 || f «-i|| + ( cfc/lT+1 + cfc3 ) + ckhT+1 + ck3 

VI - ck / VI c/c / 



< • • • 

t=i 

< ciif 1 !! + c/i T+1 + c fc 2 , 

iv 



Y^f| j -» e c as iV -> 00 and thus ( J is bounded). 
Thus, we get replacing £ n+1 , E 1 



r +l _ fi (n+l) || < c || _ || + c/l r+l + cfc 2_ 
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So we take 

(3.29) ||<T +1 || - < cp 1 -B]p\\ +ch T+1 +ck 2 . 

By (|535jl \\B^ n+1) \\ < ch T+1 + ck 2 , thus (1535)1 together with 1(535)1 . Lemma GET] for n = 0, and Lemma 
13.111 gives using the estimates of and ||£?i°' > | 

||0 n+1 || < c|| 6* 1 - S^ll + ch T+1 + ck 2 

(3.30) < cp 1 - Bf ] \\ + c||£^ 0) - B^\\ + ch T+1 + ck 2 

< cp 1 1| + ||sj 0) || + ch T+1 + ck 2 < ch T+1 + ck 2 . 



□ 



We are now ready to prove the main error estimate of this section: 
Theorem 3.13. IfO(k) = 0(h) then 

(3.31) \\U n -u(r n )\\ <ch T+1 +ck 2 , 0<n<N. 

Proof. Obviously, using Lemma T3. 121 and the fact that 9° = 0, it follows that 

\\U n - u(r n )\\ < \\e n \\ + ch T+1 < ch T+1 + ck 2 . 

4. Global Elliptic Regularity 



□ 



In this section, we present a general Global Elliptic Regularity Theorem for complex elliptic operators 
with mixed Dirichlet- Robin boundary conditions, in rectangles of M 2 . Our proof follows that of [18] which 
deals with the Dirichlet problem for real operators. In our approach, the main idea is that if the trace 
terms in the weak formulation of the problem vanish due to the boundary conditions, for suitably chosen 
test functions, then a Global Elliptic Regularity result is proved in Theorem 14.11 Note that the Robin 
condition in this Theorem does not involve any zero order term, while the first order terms are related 
to the coefficients of the boundary problem so that indeed in the weak formulation, after integration by 
parts, the trace integrals vanish. Our result is established by using the fact that the closure of a rectangle 
can be covered by using a finite union of half-balls together with an open smooth domain in the interior. 
We then apply an exponential transformation and extent our result, in Theorem l4.31 where an arbitrary 
zero order term is introduced at the Robin condition of Theorem 14.11 



Theorem 4.1. Let W = (0, 1) x (^1,^2) be a rectangular domain in cartesian coordinates. We consider 
the following boundary value problem: We seek a complex-valued function u such that 

Au zz + Bu z g + Cugg + Du z + Eu s + Fu = f in W, 

u(0,9) = 0, 

(4 1) 

V ; u{z, 6-C) = u(z, 6 2 ) = 0, 

a(9)u z + b{9)ug = at z = 1, 

where A,B,C G C^W), D,E,F G L°°{W), f G L 2 {W) and a.b : [6 1 ,6 2 ] -> C*. We also assume that 
A,B,C take imaginary values and 4-> y are always positive (or always negative). Moreover, we assume 
that 

(4.2) \AC\>^, forany(z,9)eW, 

If u £ H 1 (W) is a weak solution of (|4.1ll then the following elliptic regularity estimate holds 
(4.4) ueH 2 (W) and \\u\\ H 2 (W) < c\\f\\ L 2 (w} . 
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Proof. We consider the rectangle W. Obviously its boundary is the union of four linear segments and 
we write dW = Uf =1 dWi (cf. Figure 3). Let U z = B°{k l ,r i ) nW, be a half-ball in M 2 in W laying at 
dW of range r% and of diameter in dWi- We define its boundary by dUi := dUit U dUi C , where dUih is 
the diameter such that dUih Q dWi, and dUi C is the semicircle of range r» such that Ui C W, we also 
consider Vj = B°{ki,n/2) f] W, the half-ball being of the same center fcj as i/j and of range r</2 (cf. 
Figure 2). Obviously, cW is compact, thus dW may be covered by using a finite union of sets of the 
form Vj, while the same union together with a suitably chosen smooth domain in W covers W. By [18] 
an interior regularity estimate holds. Our aim is to prove the regularity estimate 

( 4 - 5 ) IM|jJ 3 (V<) < c||/||i,3(w 4 ), i = l,...,4. 

Interior regularity combined with the estimate (|4.5I) gives the desired result (|4.4I) (cf. [TH), pg. 322). 




FIGURE 2. Half-balls, curved boundary, horizontal boundary. 
We consider fa £ H 1 ^) and let u be the weak solution of (14. ip . If (u, v)u i :— J u uvds then we have 

( ( B \ / B \ i 

(f,<t>i)Ui =~ (Au z ,d z fa) Ui ~ u Zl d e <t>ij u + {^ Ue ^ dz( t )i ) u } 

(4.6) - (Cue,defa)u, + (Du z , fa) Ui + (Eu B , fa) Ul + (Fu, (fo)^ 



+ 



" ^' 2") +Ue (f ' C 



rjtds, 



where £>, -E are the resulting terms after integration by parts, and rfi is the outward unit normal to dUi. 
We let Cli(u, fa) := J au \u z (A, -§) + u e("§: C)]<f>iVids, and define the vector .fQ := [u z (A, ^)+ue(^,C)]fa; 
here (•, •) denotes a vector of 1R 2 . Then for dUi — dUih U dUi C it holds that £li(u, fa) — J du . h Kilfids + 
J gu Kilfids. Using the boundary conditions of u £ H 1 (W) wc obtain 



Sli(u, fa) = - / Au z fads+ / Kilfids, il 2 {u 7 fa) = / K 2 rfeds, 

JdU lh JdU lc JdU 2 c 



(4.7) Q 3 (u,fa) = - Cue fads +/ if 3 ^s, 

Q4,(u, fa) = / Cugfa\ds + / K^rjlds. 

JdUih JdU 4 c 

Our aim now is to find test functions fa such that in the weak formulation the trace terms vanish. 

Assumption 1. We assume that there exist functions fa that satisfy the following requirements: 

• The test functions are smooth and along the curved boundary Ui C oiUi vanish: fa £ H (Ui), and 
fa = Q\dU ic , i = l,...,4. 

• For i = 1,3,4, the test functions vanish also along the horizontal boundary U L h oiUi. fa = at 
z — 0, fa is arbitrary, fa = at = #i, fa = at 6 = 62- 
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Under this assumption, the sum of trace integrals in the weak formulation equals zero because Q,i(u, 4>i) 
for any i = 1, . . . , 4. The weak formulation (|4.6|) for B(u, 4>i)iAi ■— (/> 4>i)Ui becomes 

( B \ / B \ i 

. , B(u,(t>i)ui =~ {Au z ,d z fa) Ui - \(— u z ,d e 4>i) + ( — u B ,d z fa) \ 

(4.8) ^ ' w « ^ ^ ' Ui> 

- {Cu e ide4>i)ui + (Du z ,^i)ui + (Eue,4>i)ui + (Fu,<j>i) Ui - 



dW 4 



dWi 



dW 2 



9W 3 



Iz 



FIGURE 3. The rectangular domain W. 

The next step is to define, properly, for any i = !,■■■ ,4, test functions fa satisfying this assumption. 
We define the following general cut-off function ([15]) 

{0 in M 2 — B(l,r), 
1 in B(Z,r/2), 
< J < 1 elsewhere (with J = near cW c ). 

Here W := B°(l,r) n W is a half-ball in M 2 of radius r and of center [ such that dU h C iW. Let V be 
the half-ball in K 2 of center I and of range r/2 with diameter in dUh- Obviously the cut off function J 
in V equals 1, and near dU c is 0. Let u be a function in H 1 (W) that satisfies the boundary conditions of 
problem (|4.1I) . we define the function ([18]) 

(4.10) 5 := -D- h (J 2 D h u), with D fc fi(a) := gfe±MzgW , x £ U, 

where h is a positive number and e is a unitary vector (direction) in K 2 parallel to the diameter of the 
half-ball U. 

In this way for every boundary line (i = 1, ••■ ,4) of the rectangular domain W we define a cut-off 
function Ji and denote by e$ the unitary direction of the specific boundary line dWi- We then prove first 
that frj defined by these Ji in (14. 10)) for the directions are test functions that satisfy the Assumption 
1, and in the sequel we set fa :=Vi- 

More specifically, for every i = 1, ... ,4 we consider Ui = B°(ki, r,) PI W, H = 5°(^,j)n>V, fej, r, 
such that Ui C W, C 3Wi and define the cut-off function 

fj 4 = 0in M?-B(ki,n), 
Ji ■= I Ji = I in B(k h f ), 

I < Ji < 1 elsewhere (with = near dUi C ). 

Let i be a function in i/ 1 (W) that satisfies the boundary conditions of problem (|4.ip . we define as 
previously the function 

(4.11) ^ := -D7 h (J?D?u), with D*u(x) := * e Z4 
By [T5J, for any x £ Ui, the following identity holds 

(4.12) Vi(x) = — ^(Jf(x — he l )[ii(x) — u(x — he,)} — Jf(x)[u(x + hei) - u(x)]). 

Using the boundary conditions of the elliptic problem and the identity (I4.12[) . we will prove that Vi 
satisfy Assumption 1 for any i = 1, . . . , 4. 
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If i = 1 , then obviously v\ is in ff x (Wi). We notice that if x is in dlA\ c then J\(x) — and for h 
small enough J±(x — he\) = so by (14. 12)) v\(x) — {)\dU\ c . Along the boundary line dU\h holds that 
z = and e x = (0, 1). If x = (0,0) then u(x) = and u{x±he 1 ) = u(O,0±h) = 0, thus by KWi follows 
that Si (0,0) = 0. 

If i = 2 , then WaC^) e H 1 (U 2 ), and e2 = (0, 1). If x € 9^2 C then for /i small J 2 (x) = J 2 (x — he 2 ) = 0, 
thus by (|4~12l) v 2 (x) = 0\dU 2c . 

If z = 3 , then v%{x) G H 1 ^^) and e3 = (1,0), for h small. If x £ U. 3c then J 3 (x) — J 3 (x — he^) — thus 
v a (x) = 0\dU 3c . For x = (z,0i) then u{x) = u(z,0i) = and u{x ± he 3 ) = u{z ± h,0i) = 0. By flUE} 
follows that €3(2;, 6*i ) = 0. 

If z = A , then £4(2;) € H 1 ^^ and e4 = (1,0), if x is in 9^4 C then Ji(x) = and for h small enough 
J 4 (x - he 4 ) = 0, thus by (|4~T2|) f5 4 (af) = 0|5W 4c . If x = {z,6 2 ) then it (at) = u{z,0 2 ) = and 
it(a; ± hei) = u(z ± h, 2 ) = 0, thus va(z, 9 2 ) = 0. 

Therefore, in all cases Assumption 1 holds and the trace terms vanish from the weak formulation of 
the elliptic problem. If we set u := u, where u is the weak solution of the elliptic problem satisfying the 
boundary conditions, then it can be easily proved (for details see [B] and [T5]) by use of ellipticity, the 
weak formulation and the boundary conditions at z — 0, 9 = 0±, 9 = 9 2 , that for every half-ball Vj it 
holds 

( 4 -!3) ||u|| H 2 (v .) < c[||/|| i2(w .) + \\u\\ H i {Ui) ]. 

Finite summation of (|4.13p over any Vi (of type i = 1, • • • , 4) and the interior regularity give (|18j) 
(4-14) \\u\\ H 2 (w) < c[||/|| L 2 (w) + ||u|| ff i (W )]. 

Combining (|4.14j) with ellipticity we obtain the elliptic regularity result 
(4-15) Nlff 2 (W) < c||/||l2(w)- 

□ 

Remark 4.2. We note that an analogous result is also valid if in the assumptions of Theorem 14. 1[ the 
homogeneous condition at z = 1 is replaced by the non-homogeneous condition a(9)u z + b(9)ug = g at 
z = 1, for any g 6 H 2 (9W_r), where <9Wr = {1} x (9i,0 2 ). In this case, in the weak formulation the 
trace integral term containing g is hidden due to ellipticity, leaving at the right-hand side of ()4.15|) the 
extra term c\g\i 9Wr where \g\i 9Wr := inf ||u||i, for H w := {u e H X (W) : w| z=0 = 0}. More 

2 ' 2 ' v^Hw.v\ z = i— g 

specifically, the following elliptic regularity estimate holds 

(4-16) \\u\\ H 2 {w) < c||/||i2 (V v) + c\g\i^ dWR . 

The following theorem extends Theorem 14. II in the sense that we can add at the boundary condition 
along z = 1 a zero order term multiplied by an arbitrary smooth function c{9). 

Theorem 4.3. Under the assumptions 0/ Theorem \4-l\ if the boundary condition of at Z = 1 has 

the form 

(4.17) a(9)u z + b{9)u e + c{9)u{6) = at z = 1, £ [0 U 2 ], 

with c a smooth complex function of 9, then the results of Theorem \4.1\ hold (elliptic regularity). 

Proof. We set q — q(z,9) and consider the elliptic operator of (|4.ip . we apply the transformation u := 
exp(g)iu and get the following equivalent problem 

Aw zz + Bw z e + Cwee + D w w z + E w w s + F w w = f w in W, 

10(0,0) = 0, 

(4A8) w(z,0 1 ) = w(z,0 2 ) = O, 

a(0)w z +b(9)w e + c w (0)w{0) =0 at 2 = 1, 
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Table 1. Errors E{r) and spatial convergence rate for fc _1 = 400 





r = 0.1 


r = 0.5 


r = 1.0 




E{r) Rate 


E(r) Rate 


E(r) Rate 


10 


3.5162(-2) 


4.9653(-2) 


7.8266(-2) 


20 


7.5323(-3) 2.22 


1.0734(-2) 2.21 


1.6921(-2) 2.21 


40 


1.7219(-3) 2.13 


2.4518(-3) 2.13 


3.8920(-3) 2.12 


80 


4.0438(-4) 2.09 


5.7998(-4) 2.08 


9.2042(-4) 2.08 


160 


9.7655(-5) 2.05 


1.4100(-4) 2.04 


2.2381(-4) 2.04 



where D w = 2Aq z + Bq g + D,E W = Bq z + 2Cq e + E, f w = exp(-g)/, F w = F + A(q zz + q 2 z ) + B(q z g + 
Qzqe) + C(qg0 + qg) + Dq z +Eqe, and c w (9) — a{9)q z + b(8)qg + c(9). We chose q(z, 6) such that c w {9) = 
or equivalently 

(4.19) a(%,(l,0)+6(%fl(l,0) + c(0) = O for any 6 [6 x ,e % ]. 

The relation (|4"7P9l can be achieved as § = ^ i s rea l, for smooth and a(9),b(9) in C, [22 . Thus 
by (|4.18p and (|4. 19|) the problem is of the form covered by Theorem 14. 1[ and consequently 

w e H 2 (W) and ||w|| H 2 (w) < c\\f w \\ L 2 (w) . 

Obviously u — exp(q)w; therefore, u G H 2 (W) and ||u||_f/2(vv) < c II/IIl 2 (w)- ^ 

Remark 4.4. By using Remark 14.21 under the assumptions of Theorem 14.31 and if we impose the non- 
homogeneous condition a(9)u z + b(9)ug + c{9)u(9) = g at z = 1, for g € H^(dWn) in place of the 
homogeneous one, estimate (|4.16p follows (the proof is the same as in Theorem 14. 3|) . 

Remark 4.5. Theorem 14.11 and 14.31 or the results of Remarks 14.21 14.41 can be applied to cylindrical coor- 
dinates for r fixed when W = {(z,r,9) G K 3 }, by use of the change of variables u(z, 6) — u(z,9) with 
:= llif) = c °^' then the equivalent problem in cartesian coordinates is defined in a rectangular domain 
and satisfies the assumptions of Theorems 14.11 and 14.31 or those of Remarks 14.21 14.41 

5. Numerical experiments 

In this section we report on the outcome of some numerical experiments performed with the fully 
discrete scheme (|3 .4[) to solve the initial- and boundary- value problem (|1.6[) . In the notation established 
in Section 1, cf. (|1.6|) . we took T> ~ (0, l) 2 , r m in = 0, r max = 1, b = 0, = 1, D the identity matrix, 
A = (0, 1) and right-hand side F so that the exact solution is 

(5.1) u(r, y, 9) = e 2r y(e~ y - 1)8(1 - 6) 3 . 

Our first set of experiments concerns the experimental verification of the convergence rate of the 
scheme in the spatial variable. The measure of the error was the E{r) — \\u — XJ\ for r = nk, n = 1, 2, . . ., 
whereas for other values of r E was defined by linear interpolation. To determine experimentally the 
spatial order of convergence the approximate solution was computed for < r < 1 using a rectangular 
partition of D using N — hr 1 ranging from 20 to 160. The finite element space Sh consisted of piecewise 
polynomial functions of degree one. For these runs, very small r-steps were taken to ensure that the 
error due to the discretization in time-like variable r is negligible. The observed error was recorded at 
r = 0.1,0.5 and 1. As usual, the convergence rate corresponding to two different runs with mesh sizes 
hi, h>2 and corresponding errors E\ and E-x is defined to be \og(E\ / E2) / \og(h\ / hx) ■ The results are shown 
in Table [T] It is evident that the convergence rate of the spatial component of the error is indeed two. 

The determination of the accuracy in the time-like variable r is more delicate. We took = 20 
and computed the solution of our problem up to r = 1 for various values of k. For this fixed value of h 
we made a reference calculation with a small value of k = k TC f — /i/30. The corresponding approximate 
solution, denoted by Uh,rcf differs from the exact solution by a factor which is almost entirely due to the 
spatial discretization. We then define a modified measure of the error E* (r) as above but with the exact 
solution replaced by the reference solution Uh,ref- The results are shown in Table [2] 
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Table 2. Errors E{r) and r-convergence rate for h 1 = 20 



fc- 1 


E(r) E*(r) Rate 


144 


3.8104(-1) 3.9217(-1) 


192 


7.1839(-1) 1.7832(-1) 2.74 


240 


1.1771(-2) 1.0652(-1) 2.31 


288 


1.7638(-2) 7.1442(-2) 2.19 


600 


8.9952(-3) 
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